# ------------------------------------------------------
# ' run falsification tests
# ' random assignment of different grades’ entire peer to each person
# ------------------------------------------------------
load_packages = c('rio','here','dplyr','data.table','fst',
	'drake','moments','igraph','Matrix','car','lfe','expss','survey',
	'RStata','stargazer','lfe')

invisible(lapply(load_packages, library, character.only = TRUE))
library(parallel)
library(ggplot2)
library(hrbrthemes)
library(ggpubr)

# ------ specify your directory here
data_path = ""
project_path = ""

setwd(data_path)

# ------ run some codes to produce processed data sets 
source(file.path(project_path,'rdata_util.R'))

# ---- depression indicators
all_dep_indicators = c("s60i","s60j","s60k","s60l","s60m","s60n")
    #s60i : did you have trouble eating, or a poor appetite? $$$
    #s60j : did you have trouble falling asleep or staying asleep?  $$$
    #s60k : did you feel depressed or blue?  $$$
    #s60l : did you have trouble relaxing? 
    #s60m : were you moody? 
    #s60n : did you cry a lot?  $$$
    # [[ 0 : never, 1:rarely, 2:occasionally, 3:often, 4:everyday, 9:multiple response ]]

control.varlist = c('female','white','black','hispanic','other',
	'immig_1st','immig_2nd','immig_3rd','family_two','family_one','family_other','pa_educ')
    ind_controls = c('female','black','hispanic','other','immig_1st','immig_2nd',
      'family_step','family_one','family_other','pvt','pa_educ','assistance','sibsize')
    peer_controls = c(paste0(c('female','black','hispanic','other','immig_1st','immig_2nd',
      'family_one','family_other','pa_educ'),'.mean'))

processed_data_wave1 = data.table(import(file.path(data_path,'processed','v_wave1.dta')))
processed_data_wave2 = data.table(import(file.path(data_path,'processed','v_wave2.dta')))
processed_data_inschool = data.table(import(file.path(data_path,'processed','v_inschool.dta')))

raw_data_inschool = impute_aid_inschool(data.table(import(file.path(data_path,'rawdata','inschool.dta'))))
weight_wave1 = data.table(import(file.path(data_path, 'rawdata','weight_wave1.dta')))
weight_wave2 = data.table(import(file.path(data_path, 'rawdata','weight_wave2.dta')))

# additional cleaning 
cleaned_data_inschool = process_inschool(processed_data_inschool)

# process data for edge creation
processed_data_depress = process_depress(raw_data_inschool, dep_indicators=all_dep_indicators)

# --- create randomly generated 1000 regression files
edge_measure_random_entire_peer(dep_data = processed_data_depress, n_core=4)

# --- run fixed effects models
dir.create(file.path(data_path,'processed','simulation_all','grade'))

model_fe_grade  = mclapply(1:1000, function(i, type='grade') {
    random_peer_depress_grade = read_fst(file.path(data_path,'processed','simulation_all',type, paste0('sim_',type,i,'.fst')),as.data.table=TRUE)
    message('now running simulated grade ',i)
    # extract only the relevant variables
    random_peer_depress_grade = random_peer_depress_grade[, -c('scid','grade','scid_grade','depress')]
    coef = run_FE_random_peer(random_peer_depress_grade , group_fixed_effect='scid')
    coef = data.table(coef_no=1:nrow(coef),coef, run=i)
    return(coef)
},mc.cores = 6)

model_fe_grade = rbindlist(model_fe_grade)
write_fst(model_fe_grade, file.path(data_path, 'processed','simulation_all','sim_grade_model_coef.fst'), 100)
